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Abstract 


In  order  to  numerically  calculate  the  strain  energy  release  rate  (SERR)  using  finite  element 
models,  an  implementation  of  the  virtual  crack  extension  technique  was  developed  within  a 
Matlab  script.  From  the  nodal  data  read  in  from  the  ascii  output  files  produced  by  LS-Dyna,  the 
SERR  was  calculated  from  contours  of  elements  which  were  positioned  around  the  crack  tip.  In 
order  to  validate  the  implementation,  a  mesh  sensitivity  study  of  a  center  cracked  panel  was 
carried  out.  This  study  concluded  that  the  SERR  implementation  was  relatively  independent  of 
mesh  size  and  shape. 


Resume 


Nous  avons  integre  la  technique  d’extension  virtuelle  des  fissures  dans  un  script  Matlab,  afm  de 
calculer  le  taux  de  dissipation  de  Tenergie  de  deformation,  a  l’aide  d’un  modele  d’elements  finis. 
A  partir  des  donnees  nodales  tirees  des  fichiers  de  sortie  ASCII  du  logiciel  LS-Dyna,  nous  avons 
calcule  le  taux  de  dissipation  d’energie  de  deformation  pour  les  contours  des  elements  disposes 
autour  de  l’extremite  de  la  fissure.  Cette  integration  a  ete  validee  par  une  etude  de  sensibilite  pour 
le  maillage  d’un  panneau  fissure  au  centre.  La  presente  etude  a  demontre  que  notre  integration  du 
calcul  de  dissipation  de  Tenergie  de  deformation  etait  relativement  independante  de  la  taille  et  de 
la  forme  des  mailles. 
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Executive  summary 


Virtual  Crack  Extension  Technique 

K.  Shahin;  T.S.  Koko;  DRDC  Atlantic  CR  2008-274;  Defence  R&D  Canada  - 
Atlantic;  December  2008. 

Introduction  or  background:  This  report  describes  the  derivation  and  development  of  a  Matlab 
program  that  computes  the  strain  energy  release  rate  (SERR)  in  two  dimensional  finite  element 
models.  The  program  was  developed  to  post-process  displacement  results  from  LS-Dyna  using 
the  continuum  mechanics  formulation  of  the  virtual  crack  extension  technique  (VCET).  Brief 
descriptions  of  the  VCET  and  other  methods  of  extracting  SERR  results  are  summarized  as  well. 

Results:  The  SERR  for  a  center  crack  panel  are  computed  for  various  FEM  mesh  designs  and 
compared  with  analytically  derived  formulas.  The  examples  illustrate  that  the  derived  SERR 
formulations  is  correctly  implemented  for  plane  stress  models. 

Significance:  These  SERR  post-processing  routines  allow  fracture  mechanics  variables  to  be 
computed  directly  from  numerical  analyses.  By  comparing  the  computed  SERR  to  those  acquired 
experimentally,  the  onset  of  fracture  can  be  estimated  directly  from  the  simulation. 

Future  plans:  Future  plans  include  the  extension  of  the  SERR  technique  to  include  elastic-plastic 
material  behaviour  and  solid  elements.  This  will  allow  a  variety  of  geometries  and  materials  to  be 
modelled. 
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Sommaire 


La  technique  d’extension  virtuelle  des  fissures 

K.  Shahin;  T.S.  Koko;  DRDC  Atlantic  CR  2008-274;  R  &  D  pour  la  defense 
Canada  -  Atlantique;  Decembre  2008. 

Introduction  ou  contexte  :  Dans  ce  rapport,  nous  decrivons  la  production  d’un  script  Matlab 
pour  le  calcul  du  taux  de  dissipation  de  l’energie  de  deformation  pour  les  modeles 
bidimensionnels  d’elements  finis.  Ce  programme  applique  la  technique  d’extension  virtuelle  des 
fissures  avec  le  formalisme  de  la  mecanique  des  milieux  continus,  sur  les  resultats  de 
deplacements  calcules  avec  le  logiciel  LS-Dyna.  Nous  decrivons  brievement  la  technique 
d’extension  virtuelle  des  fissures  et  d’autres  methodes  d’extraction  des  resultats  du  taux  de 
dissipation  de  l’energie  de  deformation. 

Resultats  :  Nous  avons  calcule  le  taux  de  dissipation  de  l’energie  de  deformation  pour  differents 
schemas  du  maillage  du  modele  d’elements  finis  et  les  avons  compares  avec  les  formules  derivees 
analytiquement.  Ces  exemples  montrent  que  les  formulations  derivees  de  la  dissipation  de 
l’energie  de  deformation  sont  correctement  integrees  aux  modeles  de  contraintes  dans  les 
structures  planes. 

Importance  :  Ces  routines  post-traitement  de  derivation  de  la  dissipation  de  l’energie  de 
deformation  permettent  de  calculer  directement  les  variables  de  fracture  mecanique  a  partir  des 
analyses  numeriques.  On  peut  estimer  directement  le  point  de  declenchement  de  fracture,  en 
comparant  les  taux  de  dissipation  de  l’energie  de  deformation,  aux  valeurs  experimentales. 

Perspectives  :  Les  recherches  futures  porteront  sur  le  developpement  de  la  technique  de  calcul  du 
taux  de  dissipation  de  l’energie  de  deformation  pour  y  inclure  le  comportement  des  materiaux 
elastiques-plastiques  et  les  elements  solides,  ce  qui  permettra  de  modeliser  diverses  geometries  et 
substances. 
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1  Finite  Element  Formulation 


The  total  strain  energy  release  rate  (SERR)  using  the  virtual  crack  extension  method  is  given  by 
Equation  1  for  the  case  of  a  cracked  structure  with  no  body  force  component  and  in  the  absence 
traction  forces  at  the  crack  face. 


Gt 
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A  Jv 
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where, 

Gt  =  total  strain  energy  release  rate, 
Ac  =  area  of  crack  extension, 

Gij  =  stress  tensor, 

Uj  =  nodal  displacement, 

W  =  strain  energy  density, 

Sik  =  Kronecker’s  delta,  and 
Ax  =  nodal  virtual  extension, 


Using  Gauss  Quadrature  integration  method,  the  SERR  is  obtained  using  the  following  expression 
to  approximate  the  volume  integral  given  by  Equation  1  for  two-dimensional  plane  stress  or  plane 
strain  conditions: 
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where, 

Aa  =  virtual  crack  extension, 

Ne  =  number  of  elements  in  the  model, 
n,  m  =  order  of  Gauss  Quadrature  integration, 

[/]  =  2x2  identity  matrix, 

{U0}  =  nodal  displacement  vector, 

{AX0}  =  virtual  extension  of  element  nodes, 
[P]  =  derivatives  of  shape  functions, 

[J]  =  element  Jacobian,  and 
a  =  Gauss  quadrature  weight  factors. 
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See  deLorenzi  (1985)  for  more  detailed  description  of  the  expressions  shown  above  and  their 
derivations. 


Instead  of  relying  on  stress  and  strain  outputs  from  LS-Dyna,  to  both  reduce  the  amount  of 
manual  pre-processing  required  and  to  maintain  consistency  in  the  formulation,  the  required  stress 
and  strains  are  derived  from  LS-Dyna  displacement  results  as  given  below 


a 


dU 


e\  =  < 


dX 

dU 


(1.1) 

(2,2) 


dX 


(3) 


Therefore,  nodal  displacement  results  are  the  only  required  output  from  LS-Dyna,  and  only  for 
the  elements  involved  in  the  SERR  calculations.  Furthermore,  the  SERR  contribution  of  elements 
in  which  all  nodes  undergo  the  same  (or  no)  crack  extension  is  exactly  zero.  Gauss  quadrature 
integration  is  performed  using  one -point  rule,  in  consistency  with  the  LS-Dyna  reduced 
integration  of  the  element  stiffness  matrix. 
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2  Procedure: 


To  use  the  vcet_martec.m  program,  the  following  procedure  is  recommended. 

Define  the  elements  surrounding  the  crack  under  a  unique  part  ID  (component)  separate  from  the 
rest  of  elements  in  the  model.  These  elements  will  be  used  in  the  calculation  of  the  SERR  in  the 
vcet  martec.m  program. 

For  convenience,  re-number  the  nodes  and  elements  in  this  component  to  start  at  1,  and  the  rest  of 
the  elements  and  nodes  in  the  model  to  start  a  much  larger  number.  This  is  not  strictly  necessary, 
however,  it  greatly  improves  the  experience  of  specifying  the  nodal  shifting  to  be  specified  by  the 
user  and  the  process  of  obtaining  element  nodal  coordinates  and  displacements. 

The  direct  input  to  the  program  is  done  through  the  text  file  input_parameters.dat,  which 
contains  the  following  comma  separated  input  parameters: 

filel  .csv,  file2.csv,  file3.csv,  file4.csv,  psID,  symmetrylD,  cwID,  modulus,  pr,  deltax 

where 

filel. csv  :  contains  the  node  coordinates  input  (character  string) 
file2.csv  :  contains  the  element  nodal  connectivity  input  (character  string) 
file3.csv  :  contains  the  nodal  displacement  solution  (character  string) 
file4.csv  :  contains  the  nodes  affected  by  virtual  crack  extension  (character  string) 

psID:  Element  formulation  option:  (integer) 

EQ.  12:  Plane  stress  (default) 

EQ.  13:  Plane  strain 

symmetrylD:  Model  symmetry  option  (integer) 

EQ.  0:  Not  symmetric  (default) 

EQ.  1 :  Crack  lies  on  plane  of  symmetry. 

cwID  :  Element  node  numbering  order  (integer) 

EQ.  0:  Counter  clockwise  (default) 

EQ.  1 :  Clockwise. 

modulus  :  Elements  Young’s  modulus  (float) 

pr  :  Poisson  ratio  (float) 

delta  x  :  Virtual  crack  extension  (float) 

Example  of  input  line  (used  in  verification  example  1): 

c_nodes.csv,  c_elements.csv,  c_disp.csv,  c_paths.csv,  12,  1,  1,  207000,  0.3,  0.001 
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Note:  character-type  (file  name)  input  is  case  sensitive. 

As  mentioned  above,  the  vcet  martec.m  program  requires  four  input  files.  Below,  a  brief 
description  is  given  on  the  format  and  input  parameters  of  the  four  files. 

1.  filel  .csv:  nodal  coordinates  input 

From  the  LS-Dyna  input  file,  copy  the  nodal  coordinates  of  the  nodes  attached  to  the 
crack-elements  component  into  Excel  and  remove  the  column  corresponding  to  the  z- 
coordinate.  Save  this  information  in  a  comma  separated  “.csv”  format. 

Input  format:  Node  number,  x-coordinate,  y-coordinate 

2.  file2  .csv:  element-node  connectivity  input 

Obtain  a  list  of  elements  nodal  connectivity  from  the  LS-Dyna  input  file.  Copy  this 
information  into  Excel  and  remove  the  column  corresponding  to  the  Part  ID  Number 
(second  column  in  LS-Dyna  input).  Save  this  information  in  comma  separated  “.csv” 
format. 

Input  format:  Element  number,  node  1,  node  2,  node  3,  node  4 

3.  file3.csv:  node  displacement  solution 

Copy  the  nodal  displacement  solution  obtained  from  LS-Dyna  “nodeout”  ASCII  file  into 
Excel,  and  save  only  the  columns  corresponding  to  the  node  numbers,  and  displacements 
in  the  x-  and  y-  directions. 

Input  format:  node  number,  x-displacement,  y-displacement 

4.  file4  .csv:  nodes  affected  by  virtual  crack  extension 

Specify  the  nodes  affected  by  the  virtual  crack  extension.  Obtain  the  node  numbers  from 
LS-Dyna  input  file,  list  them  in  Excel  such  that  each  column  contains  the  nodes  to  be 
shifted  in  each  path,  and  save  in  comma  separated  “.csv”  format.  It  is  assumed  that  the 
nodes  undergoing  a  shift  for  a  given  path  are  also  shifted  in  all  subsequent  paths,  and 
need  not  be  specified  more  than  once.  This  reduces  the  amount  of  user  input  required  to 
specify  the  node  shifts.  To  reduce  the  programming  effort,  the  number  of  entries 
must  be  the  same  in  all  columns  (paths)  in  the  “.csv”  input  file,  therefore,  enter  0  in 
place  of  blank  node  numbers. 

Input  format:  comma  separated  matrix  of  node  numbers  [Ny],  ith  node  in  j  path. 

The  user  should  exercise  caution  when  saving  the  above  files,  especially  the  nodal-displacements 
comma  separated  “.csv”  file  to  ensure  a  sufficient  number  of  significant  figures  are  saved. 
During  the  development  phase,  it  was  noticed  that  Excel  often  displayed  the  displacement  results 
in  scientific  notation  using  only  two  decimal  places.  If  left  unchanged,  only  the  displayed 
significant  figures  would  be  saved  into  the  “.csv”  file,  and  the  lost  significant  figures  cannot  be 
retrieved. 

Initially,  the  program  was  developed  using  counter-clockwise  (CCW)  node  numbering,  consistent 
with  most  finite  element  formulations.  However,  positive  crack  opening  displacements  could 
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only  be  obtained  by  reversing  the  element  normals  such  that  nodes  are  numbered  in  a  clockwise 
(CW)  fashion.  Therefore,  CW  node  numbering  is  used  in  this  program.  It  is  interesting  to  note 
that  in  the  absence  of  contact  definitions  to  prevent  negative  crack  opening  displacements,  using 
CCW  node  numbering  and  negative  displacements  results  in  positive  and  accurate  estimates  of 
the  SERR,  which  highlights  the  importance  of  first  checking  the  deformed  shape  to  ensure 
physically  admissible  deformations. 
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3  Brief  Overview  of  Alternative  Techniques 


The  method  developed  by  deLorenzi  to  compute  SIF  uses  a  continuum  mechanics  approach  to 
apply  the  VCET.  Under  elastic  conditions,  the  approach  is  analytically  equivalent  to  the  J- 
integral  (using  Green's  theorem  to  convert  the  path  integration  to  an  equivalent  area  integral). 

The  simplest  means  of  extracting  SIFs  from  finite  elements  results  is  through  the  so-called  crack- 
tip  opening  displacement  (CTOD)  method.  From  linear-elastic  fracture  mechanics,  the  stress 
intensity  factors  are  used  to  describe  the  stress  and  displacement  fields  around  the  crack-tip. 
Therefore,  knowing  the  crack-tip  displacement  field  from  finite  element  results,  one  can  compute 
the  stress  intensity  factors  due  to  the  applied  load.  Using  two-dimensional  four-node 
quadrilateral  element  configuration  shown  in  Figure  1,  the  stress  intensity  factors  from  this 
technique  are  given  by  Equation  4. 


K, 


K, 


IjUy/lx  vb  -va 

4k  *- 

ub  —ua 

4k  k 


(4) 


where, 

u  and  v  =  axial  and  lateral  displacements,  respectively,  of  the  crack-tip  nodes  shown  in 
Figure  1, 

Le  =  length  of  the  element, 
p  =  material  shear  modulus,  and 
k  =  material  parameter  given  by 


3-4v 
k=3-v 
l  +  v 


plane  stress 
plane  strain 


It  can  be  shown  that  strain  energy  release  rate  is  related  to  the  stress  intensity  factor  through 


where,  E*  is  the  effective  elastic  modulus  given  by; 


E  in  plane  stress 
Er=\  E 

,  - in  plane  strain 

ll-v2 
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An  analogous  procedure  is  the  modified  crack  closure  integral  (MCCI)  developed  by  Rybicki  and 
Kanninen  (1977),  which  has  proven  very  popular  in  the  literature.  However,  it  is  strictly 
applicable  to  linear  elastic  conditions.  In  the  MCCI,  the  SERR  is  estimated  by  the  energy 
required  to  close  the  crack  by  a  small  amount  Aa.  In  validating  the  results  of  vcet  martec.m 
results,  comparisons  are  made  to  MCCI  results.  According  to  the  MCCI  technique,  the  SERR  in 
modes  I  and  II  are  given  by 


Gj  = 


Gu  = 


24 

K_ 

2  L, 


(Vb~Vc) 

iub-uc) 


where, 

Ff  and  Fya  =  crack-tip  nodal  sliding  and  opening  forces,  respectively,  and 
u  and  v  =  sliding  and  opening  displacements,  respectively,  of  the  crack-tip  nodes  shown 
in  Figure  1. 


Figure  1.  Finite  element  mesh  around  the  crack-tip  using  four-node  quadrilateral  elements  for 

CTOD  and  MCCI  calculations. 

The  main  advantage  of  the  CTOD  and  MCCI  over  the  VCET  method  lies  in  the  fact  the  strain 
energy  release  rate  components  are  determined  in  each  mode  separately,  whereas  the  VCET, 
similar  to  the  J-integral,  only  gives  estimates  of  the  total  SERR.  To  separate  the  contributions  of 
modes  I  and  II,  Bui  (1983)  proposed  a  procedure  to  separate  modes  I  and  II  contributions  to  the  J- 
integral  (which  can  also  be  applied  to  the  VCET)  by  separating  the  displacement,  stress  and  strain 
fields  into  a  symmetric  (mode  I)  and  an  anti-symmetric  (mode  II)  components. 

Typically,  the  VCET  and  J-integral  methods  are  more  accurate  than  the  MCCI,  and  all  are  far 
more  accurate  than  the  CTOD  method.  In  fact,  the  CTOD  can  only  produce  acceptable  results 
using  singular  eight-node  elements.  However,  the  computational  effort  required  for  programming 
the  VCET  or  J-integral  methods  for  an  arbitrary  crack  under  mixed-mode  conditions  is 
substantially  greater  than  that  required  for  the  MCCI  method. 
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3.1  Possible  Sources  of  Error  in  calcj.m  Program 


Unfortunately,  we  could  not  run  the  Matlab  code  provided  by  the  client  (calcj.m),  and  the  error 
could  not  be  readily  resolved.  However,  three  possible  sources  of  error  became  apparent  by 
inspecting  the  structure  of  the  calcj.m  code,  and  they  are: 


1. 


The  strain  energy  density  is  computed  as  half  the  product  of  the  stress  and  strain  tensors. 
This  approach  would  result  in  incorrect  estimates  of  the  strain  energy  density  since 


1  1  *„Y. 


xy /  xy 


whereas  the  tensor  product  yields  the  incorrect  value  given  below,  since  the  shear  stress 
component  is  accounted  for  twice. 
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LS-Dyna  output  provides  the  true  strain,  yxy,  which  is  twice  the  engineering  strain,  sxy. 
Therefore,  the  strain  matrix  used  in  “ calcj.m ”  must  take  this  into  account  by  modifying 
the  LS-Dyna  output  results  before  use  in  the  SERR  calculations. 


2.  The  correct  weight  function  for  one-point  Gauss  Quadrature  integration  is  2.0,  and  not 
1.0  as  used  in  “ calcj.m ”  code. 

3.  The  correspondence  between  the  element  numbers  in  “elcont”  and  those  of  “csig”  and 
“ceps”  is  not  readily  apparent.  Element  stresses  referenced  in  the  range  in  the  file 
“elcont”  number  between  5027  and  5527,  however,  only  4950  element  stress  and  strain 
tensors  were  provided  in  files  “csig”  and  “ceps”. 
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4  Verification  Examples: 


4.1  Example  1 

The  SERR  is  evaluated  for  a  centrally  cracked  1 0  mm  square  plate.  The  plate  thickness  is 
0.1  mm,  and  the  crack  length  is  0.6  mm.  Young’s  modulus  and  the  Poisson  ratio  of  the  plate  are 
207  GPa  and  0.3,  respectively.  According  to  linear  elastic  fracture  mechanics,  the  exact  SERR  is 
given  by 


where, 

a  =  half  crack  length, 
b  =  half  plate  width,  and 
s  =  remote  tensile  stress, 

The  model  is  double-symmetric,  therefore  only  a  quarter  of  the  plate  is  modeled.  The  finite 
element  mesh  used  in  this  example  is  relatively  coarse  (558  nodes)  with  a  focused  circular  15x16 
mesh  around  the  crack-tip  as  shown  in  Figure  2.  Virtual  crack  extension  is  introduced  in 
seventeen  paths,  the  first  contains  only  the  crack-tip  node,  and  each  successive  path  contains 
a  layer  of  nodes  surrounding  the  crack-tip.  A  truncated  version  of  the  LS-Dyna  input  file  is 
given  in  Annex  A,  and  the  full  LS-Dyna  input  is  included  in  the  electronic  appendix 
{circular Elements  05.  dyn ) . 

As  illustrated  in  the  truncated  version  of  the  input  file,  the  elements  surrounding  the  crack  are 
defined  under  a  separate  part  ID,  and  the  nodes  attached  to  these  elements  are  numbered  1 
through  256.  Nodes  outside  this  component  have  numbers  greater  than  1000.  Similarly,  elements 
surrounding  the  crack  are  numbered  1  through  240,  whereas  elements  outside  the  crack  region  are 
numbered  1241  and  higher. 

Electronic  copies  of  all  input  files  used  in  this  model  are  provided  in  electronic  format  in  the 
vcet_martec.zip  file.  The  input  command  line  for  this  model  is  given  by 

cnodes.csv,  celements.csv,  cdisp.csv,  c_paths.csv,  12, 1,  1,  207000,  0.3,  0.001 


In  the  first  verification  example  used,  the  average  SERR  from  all  fourteen  paths  that  include  more 
than  1  node  is  -1.44  %,  whereas  the  error  in  the  MCCI  result  is  2.85  %.  The  CTOD  cannot  be 
expected  to  yield  practically  usable  estimates  of  the  SERR  due  to  the  use  of  first  order  elements 
(which  cannot  capture  the  singular  displacement  field).  The  error  in  the  CTOD  SERR  is  247  %. 
The  results  are  summarized  in  Table  1  for  all  paths,  including  one  where  all  nodes  are  shifted 
which  results  in  zero  SERR  (serves  to  validate  the  integrity  of  the  model). 
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Lastly,  the  vectmartec.m  program  is  used  to  determine  the  SERR  in  two  models  using  triangular 
elements  (shown  in  Figure  3)  and  a  very  coarse  mesh  (shown  in  Figure  4).  The  results  from  the 
finite  element  model  are  summarized  in  Table  1,  which  illustrates  a  common  feature  of  energy 
based  methods  of  extracting  SERR  from  FE  models  in  that  an  overly  fine  mesh  is  not  necessary  to 
obtain  practical  results.  Furthermore,  the  applicability  of  the  model  is  verified  for  both 
quadrilateral  and  triangular  elements. 

4.2  Example  2 

4.2.1  Objectives 

i)  Verify  the  applicability  of  vcetmartec.m  in  finite  element  models  using  a 
rectangular  elements  of  quadrilateral  elements,  and 

ii)  Verify  the  accuracy  of  the  vcet  martec.m  program  in  instances  where  the  crack  plane 
does  not  coincide  with  a  plane  of  symmetry. 

A  10-mm  square  steel  plate  with  a  central  crack  is  analyzed  using  FS-Dyna,  and  the  strain  energy 
release  rates  (SERR)  are  determined  using  the  virtual  crack  extension  technique  (VCET)  as 
implemented  in  the  vcet  martec.m  code.  All  files  required  to  run  this  problem  are  included  in  the 
attached  rectangular_mesh.zip  file.  The  input  parameters  are  given  by: 

rectnodes.csv,  rectelements.csv,  rectdisp.csv,  rectpaths.csv,  12,  0, 1,  207000,  0.3,  0.001 

This  is  the  same  geometry  and  load  conditions  used  in  the  earlier  verification  example,  except  the 
mesh  consists  entirely  of  quadrilateral  elements  in  the  crack-tip  region  as  shown  in  Figure  5. 
Furthermore,  unlike  in  the  previous  example,  this  model  does  not  take  advantage  of  symmetry 
about  x-axis  (note,  the  input _parameters.dat  control  card  number  six  is  set  to  zero  as  given 
above)  since  the  crack-tip  is  not  defined  on  a  plane  of  symmetry  (as  shown  in  Figure  6).  The 
SERR  obtained  from  nine  paths  are  summarized  in  Table  2,  where  in  the  ninth  path  all  nodes  are 
shifted  by  an  equal  amount  to  ensure  the  resulting  SERR  is  zero  in  such  instances.  As  expected, 
least  accurate  estimates  of  SERR  are  obtained  when  the  crack  is  virtually  extended  by  the  shifting 
of  only  a  single  node. 

4.3  Element  Thickness  Effects 

The  formulation  of  the  SERR  given  in  vcet_martec.m  is  in  fact  independent  of  the  element 
thickness,  assuming  the  elements  in  the  crack  component  are  of  the  same  thickness.  It  was  noted 
that  the  calcj.m  code  incorrectly  divides  the  result  of  the  SERR  integral  by  (Aa  x  t),  when  in  fact  it 
should  only  be  divided  by  (Aa)  for  two-dimensional  plane  stress  and  plane  strain  elements.  The 
thickness  effects  are  in  fact  cancelled  out  by  performing  the  volume  integration  as  given  below 

J7  dV  J7  t  dA  JldA 

V_ _ _  A _ constant  _  a _ 

A  at  ~  A  at  thickness  ~  A  a 
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The  expressions  given  by  deLorenzi  (1985)  also  reflect  this  observation.  Specifically,  equation 
21  in  deLorenzi ’s  work  applies  to  plane  stress  and  plane  strain  elements,  regardless  of  the  element 
thickness,  as  long  as  all  elements  (in  the  integral)  share  the  same  thickness. 

However,  it  was  noted  that  LS-Dyna  pressure  input  (using  *LOAD_SEGMENT)  only  produced 
the  correct  displacement  magnitudes  when  a  unit  thickness  is  used  in  *SECTION_SHELL 
definition.  Therefore,  all  verification  examples  were  developed  using  a  unit-thickness  plane 
stress  elements.  For  elements  with  different  thickness  values,  the  program  vcet  martec.m  would 
still  produce  the  correct  results,  given  that  the  nodal  displacement  results  from  LS-Dyna  are 
correctly  obtained.  This  observation  (most  likely)  does  not  apply  to  models  when  the  load  is 
applied  through  nodal  loading. 

For  example,  in  the  attached  rectMesh_06.dyn  model  (included  in  rectangular_mesh.zip), 
changing  the  element  thickness  under  *SECTION_SHELL  from  1.0  to  0.1  would  result  in  a  ten¬ 
fold  increase  in  displacements.  This  was  not  expected,  since  the  load  is  applied  as  uniform 
pressure  using  the  control  cards  shown  under  *LOAD_SEGMENT.  It  was  expected  that  LS- 
Dyna  would  factor  in  the  element  thickness  effects  when  the  pressure  load  is  converted  to  an 
equivalent  nodal  point  load. 
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Table  1.  SERR  Estimates  from  Example  1  FEA  Models 


SERR  (N/mm) 

Path 

Circular 

mesh 

Triangular 

mesh 

Coarse 

mesh 

1 

8.2896 

8.7604 

8.2670 

2 

9.1198 

9.3749 

10.0396 

3 

9.2698 

9.2785 

0 

4 

9.3268 

9.2442 

5 

9.3607 

9.2368 

6 

9.3823 

9.2401 

7 

9.3978 

9.2434 

8 

9.4078 

9.2459 

9 

9.4165 

0 

10 

9.4271 

11 

9.4274 

12 

9.4343 

13 

9.4378 

14 

9.4387 

15 

9.4448 

16 

0 

Average 

9.305413 

9.203007 

9.153285 

MCCI 

9.5081 

9.2827 

10.3381 

CTOD 

31.1626 

22.0386 

53.1984 

E  (GPa) 

207 

Exact 

9.24870 
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Table  2.  Summary >  of  Example  2  Results. 


Path 

SERR  (N/mm) 

1 

7.4031 

2 

8.9554 

3 

9.1995 

4 

9.2632 

5 

9.2747 

6 

9.2896 

7 

9.2867 

8 

9.2858 

9 

0 

Average  (paths  1-8):  8.99475  N/mm 

Average  (paths  2-8):  9.22129  N/mm 

Exact  =  9.24870  N/mm 
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Uniform  Pressure  =  1000  MPa 


>TVt: 

t  rt 

Plane  of  symmetry 


Figure  2.  Finite  element  model  using  circular  mesh  around  the  crack-tip  (Example  1). 
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Figure  3.  Finite  element  model  using  triangular  elements  around  the  crack-tip  (Example  1). 
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Plane  of  symmetry 


Uniform  Pressure  =  1000  MPa 


Plane  of  symmetry 


Figure  4.  Finite  element  model  using  a  coarse  mesh  (Example  1). 


16 


DRDC  Atlantic  CR  2008-274 


Uniform  Pressure  =  1000  MPa 
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Figure  5.  Finite  element  model  using  a  rectangular  mesh  around  the  crack-tip  (Example  2). 
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Figure  6.  Close-up  of  elements  surrounding  the  crack-tip  (Example  2). 


18 


DRDC  Atlantic  CR  2008-274 


References 


HD  Bui,  1983,  Journal  of  the  Mechanics  and  Physics  of  Solids,  1983,  439-448 
HG  Delorenzi,  1985,  Engineering  Fracture  Mechanics,  (21)  129-143. 

EF  Rybicki,  and  MF  Kanninen,  1977,  Engineering  Fracture  Mechanics  (9)  931-938. 


DRDC  Atlantic  CR  2008-274 


19 


This  page  intentionally  left  blank. 


20 


DRDC  Atlantic  CR  2008-274 


Annex  A  Sample  Input  Fie 


*KEYWORD 

"TITLE 

Circular  Crack  Mesh 

$ 

$  Solution  parameters  not  shown 

$ 

*NODE 

1  0.5260634364773  0.0305469888013  0.0 

$ 

$  Nodes  in  the  CrackElements  component 

$  256  nodes  in  total 

$ 

256  0.58  0.0  0.0 

1231  0.2674344206075  0.461001262414  0.0 

$ 

$  Nodes  in  the  RestOfElements  component 

$  302  nodes  in  total 

$ 

1580  2.0322314049587  2.7636363636364  0.0 

$ 

$  Material  definition 

$ 

*MAT_ELASTIC 

1  7.900E-09  207000.0  0.3 

$ 

$  Definition  of  element  parts  (components) 

$ 

TART 

$ 

$  Part  ID  2  contains  the  CrackElements  component 

$ 

2  1  1 

$ 

$  Part  ID  3  contains  the  RestOfElements  component 

$ 

3  1  1 

*  SECTIONSHELL 
1  12 

0.1  0.1  0.1  0.1 

*ELEMENT_SHELL 

2  2  194  12  26  26 

17  2  194  26  40  40 

$ 

$  Triangular  elements  at  the  crack-tip,  12  elements  in  total 
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$  Note  repeated  node  definition  for  nodes  3  and  4 

$ 

211  2  194  242  228  228 

226  2  194  256  242  242 

*ELEMENT_SHELL 

1  2  26  12  13  27 

$ 

$  Quadrilateral  elements  in  the  CrackElements  component 

$  224  elements  in  total 

$  Therefore,  total  number  of  elements  in  the  CrackElements 
$  component  is  240  elements 

$ 

240  2  243  209  216  229 

1241  3  1231  1249  1260  1261 

$ 

$  Elements  in  RestOlElements  component. 

$ 

1519  3  1580  1572  1558  1559 

$$ 

$$  Sets  Defined  In  ElyperMesh 

$$ 

*SET_NODE_LIST 

1 

194 

*SET  NODE  LIST 


2 

12 

26 

40 

54  68  82  96  110 

124 

138 

165 

166  193  195  228  242 

256 

$ 

$  Total  of  16  node  sets  for  the  paths  defining  the  crack  extension 

$ 

*SET_NODE_LIST 

16 

209  210  211  212  213  214  215  216 

217  218  219  220  221  222  223  224 

225 

*SET_NODE_LIST 

101 

$ 

$  Contains  nodes  1  to  256  (defined  in  16  sets  above)  in  one  set ...  this  set 

$  is  used  to  save  node  displacement  results  required  for  fracture  analysis. 

$ 

*DAT ABASE  HISTORY  NODE  SET 
101 

*END 
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